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CHAPTER  1 


HIGH  RESOLUTION  TWO  DIMENSIONAL  SHEAR  BAND 
COMPUTATIONS:  IMPERFECTIONS  AND  MESH  DEPENDENCE 

1.1.  INTRODUCTION 

The  calculation  of  shear  bands  and  other  material  instability  phenomena  has  become 
of  considerable  interest  because  of  its  importance  in  predicting  material  failure.  This  paper 
deals  with  two  issues  for  this  class  of  problems: 

1.  the  role  of  imperfections  in  setting  the  length  scale  and  structure  of  shear  bands 
in  two  dimensional  problems; 

2.  the  behavior  of  various  elements  types  for  shear  bands  which  are  not  aligned 
with  element  edges. 

The  primary  tool  used  in  examining  these  issues  is  a  massively  parallel  version  of  an 
explicit  program,  Belytschko,  Plaskacz,  Kennedy,  Greenwell  (1990),  and  Belytschko, 
Plaskacz,  Chiang  (1991).  By  means  of  this  program,  we  were  able  to  solve  problems 
using  up  to  64,000  elements. 

For  coarse  meshes,  numerical  results  are  strongly  mesh  dependent  for  the  shear 
band  problem.  Tvergaard  etal  (1981)  and  Needleman  and  Tvergaard  (1984)  have  used 
crossed-diagonal  meshes  of  constant  strain  triangles  and  J2  flow  theory  with  strain- 
softening  to  calculate  shear  bands.  The  crossed-diagonal  meshes  developed  by  Nagtegaal, 
Parks  and  Rice  (1974)  were  used  to  avoid  volumetric  locking.  For  these  elements  and  the 
refinement  which  was  used,  it  was  necessary  for  the  mesh  to  be  arranged  so  that  the 
diagonals  of  the  mesh  coincided  with  the  orientation  of  the  shear  band.  Otherwise, 
localization  would  not  occur.  This  behavior  of  numerical  solutions  is  often  called  mesh 
dependence. 

Several  approaches  have  been  used  to  decrease  mesh  dependence  for  shear  bands. 
Ortiz,  Leroy  and  Needleman  (1987)  used  the  loss  of  ellipticity  in  an  element  to  trigger  the 
introduction  of  special  step-function  strain  fields  associated  with  incompatible  modes  of  the 
quadrilateral  element.  Fish  and  Belytschko  (1988)  used  embedded  band-like  functions 
triggered  by  loss  of  ellipticity  to  obtain  good  resolution  with  coarse  meshes.  Pietruszczak 
and  Mroz  (1981)  first  proposed  a  similar  approach  for  triangular  elements.  However,  these 
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approaches  can  only  give  improvement  in  the  overall  response  of  a  body,  such  as  the  force- 
deflection  curve.  The  detailed  structure  of  the  strain-field  within  the  shear  band,  including 
the  maximum  strain,  cannot  be  computed  from  these  approaches. 

Belytschko,  Fish  and  Bayliss  (1990)  used  a  spectral  overlay  on  finite  elements, 
which  was  able  to  provide  a  detailed  resolution  in  the  shear  band.  They  showed  the  critical 
role  of  the  initial  imperfection  on  the  morphology  of  the  strain  field  in  the  shear  band. 
However,  the  method  entailed  a  priori  knowledge  of  the  location  of  the  shear  band. 
Recently,  an  r-adaptive  mesh  refinement  scheme  has  been  applied  to  localization  problems 
by  Ortiz  and  Quigley  (1991)  and  it  proved  quite  effective,  but  detailed  results  on  shear  band 
structure  were  not  given. 

The  critical  issue  of  length  scales  in  solutions  of  viscoplastic  problems  has 
perplexed  numerous  researchers.  Quasi-static  viscoplasticity,  in  the  absence  of  heat 
conduction,  lacks  a  length  scale,  and  in  dynamic  problems  the  length  scale  is  markedly 
smaller  than  the  dimensions  of  shear  bands.  Several  investigators,  such  as  Pan  (1983)  and 
Shawki  and  Clifton  (1989),  have  mentioned  imperfections,  but  their  critical  role  in 
determining  the  morphology  of  a  shear  band  has  not  been  studied  extensively.  Needleman 
(1988)  demonstrated  the  role  of  imperfection  size  for  stepwise  imperfections  in  one- 
dimension.  Belytschko,  Bayliss  and  Fish  (1990),  by  using  a  spectral  overlay  on  finite 
elements  which  provided  higher  resolution  at  the  shear  bands,  were  able  to  show  that  the 
length  and  shape  of  an  imperfection  plays  a  critical  role  in  determining  the  length  scale  of  a 
shear  band,  as  exhibited  by  its  width.  Subsequently,  Belytschko,  Moran  and  Kulkami 
(1991)  developed  a  closed-form  solution  for  a  class  of  one-dimensional  viscoplastic 
problems  which  clearly  demonstrated  the  dominant  role  of  the  imperfection  in  determining 
the  shear  band  structure  after  the  onset  of  material  instability.  They  showed  that  the  initial 
width  of  the  shear  band  is  completely  determined  by  the  initial  width  of  the  imperfection. 
Furthermore,  they  showed  that  step  function  imperfections,  which  are  widely  used  in  two- 
dimensional  finite  element  calculations,  are  of  questionable  value  insofar  as  the  morphology 
of  the  shear  band  is  concerned  because  these  lead  to  unstable  solutions.  Batra  and 
Liu(1990)  have  also  studied  the  role  of  imperfections  on  adiabatic  shear  bands  in  two 
dimensions. 

In  this  paper,  dynamic  formation  and  progression  of  two-dimensional  shear  bands 
isstudied  using  much  finer  meshes  than  have  previously  been  reported.  Various  sizes  of 
imperfections  are  considered  and  both  crossed-diagonal  triangular  element  meshes  and 
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quadrilateral  meshes  are  used.  It  will  be  shown  that  the  onset  of  localization  and  the 
structure  of  the  shear  band  in  the  early  stages  of  localization  is  independent  of  element  type 
if  the  mesh  is  sufficiently  refined  compared  to  the  scale  of  the  imperfection.  However,  in  a 
viscoplastic  shear  band  the  strain  is  characterized  by  an  exponential  growth  which  depends 
on  the  magnitude  of  the  initial  imperfection,  so  that  the  region  of  highest  strain  gradient 
narrows  with  time,  and  eventually  only  a  few  elements  span  this  region.  At  this  point,  the 
solution  becomes  element  and  mesh-dependent. 

It  is  also  shown  that  in  two  dimensions,  as  in  one  dimension,  the  structure,  location 
and  size  of  the  imperfection  play  a  crucial  role  in  the  formation  of  the  shear  bands.  For 
larger  imperfections,  the  initial  shear  band  is  markedly  wider,  and  in  these  large-scale 
computations  often  spans  6  to  10  elements.  An  interesting  transition  in  the  arrangement  of 
shear  bands  was  also  noted  in  the  computations:  when  the  size  of  the  imperfection  is 
reduced  in  one  example,  a  transition  to  a  set  of  parallel  shear  bands  was  noted.  Such 
transitions  have  also  been  observed  in  experiments. 

The  paper  is  organized  as  follows:  Section  1.2  gives  a  brief  outline  of  the 
computational  method.  Section  1,3  describes  the  computations  and  the  results,  and  Section 
4  gives  some  remarks  and  conclusions. 

1.2.  NUMERICAL  METHOD 


The  problems  were  solved  by  a  finite  element  method  with  a  Lagrangian  mesh.  The 
rate-of-deformation  tensor 


or 


T|  =  D  V 


(1) 


is  used  as  a  measure  of  rate-of-deformation;  here  v^  are  the  velocities,  Xj  the  Eulerian 
(spatial)  coordinates. 


The  constitutive  equation  was  expressed  in  terms  of  the  Jaumann  rate  of  Cauchy 

stress 


a  =  Cq  =  CD  V 

.  V  T 

0=0+  wo  +  ow  ^ 


(2a) 


3 


(2b) 


where  a  is  the  Jaumann  rate,  C  is  a  constitutive  tangent  matrix  which  depends  on  the 
current  state  of  stress,  and  w  is  the  spin  tensor. 

The  governing  equations  are  Eqs.  (1)  and  (2)  and  the  momentum  equation 

tJijj  +  hi  =  pvi  Qj.  jy  +  b  =  pv  (3) 

In  two  dimensions,  Eqs.  (1-3)  represent  5  equations  in  the  5  unknowns  Sjj  and  vj. 

The  finite  element  equations  are  obtained  by  approximating  the  velocities  in  each 


element  by 

v(X,t)  =  N(X)de(t)  (4) 

The  resulting  ordinary  differential  equations  are; 
momentum  equation 

Ma(t)  =  f(t)  a  =  =  d(  t)  (5) 

f(0  =  fext(0  -  fint(0  (6) 

^  Lel«xt(0  fint=  ^  Lefint(0 

*  e  (7) 

^nt(t)=|  B^adQ  (8) 

Jih 

measure  of  rate-of-deformation 

Ti(X,t)  =  B(X)de(t)  B=Z>N  (9) 

de(t)  =  Led(t)  (10) 

and  the  constitutive  equation  (2).  In  the  above 
t  =  time; 

M  =  mass  matrix; 
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d  =  nodal  displacement  matrix; 

X  =  material  coordinate; 

f(t).  fext(0.  *int(0  =  resultant,  external  and  internal  nodal  forces,  respectively; 

B(X)  =  semidiscrete  form  of  the  symmetric  part  of  the  gradient  operator; 

11  =  rate-of-deformation  (velocity  strain)  tensor; 

Le  =  Boolean  connectivity  matrix  of  element  e. 

The  SIMD  algorithm  is  described  in  Belytschko,  Plaskacz,  Kennedy  and  Greenwell 
(1990).  The  algorithm  uses  an  exchange  of  forces  at  each  time  step  and  integrates  the 
nodes  redundantly  so  that  only  one  communication  is  required  in  each  time  step,  in  contrast 
to  the  two  communications  required  in  a  standard  gather-scatter  finite  element  algorithm. 
The  exchange  algorithm  can  also  be  used  to  advantage  in  MIMD  computers.  The  essential 
idea  here  is  to  use  the  gather-scatter  operations  within  each  subdomain,  so  that  gather- 
scatter  does  not  require  communication  between  processors.  The  exchange  operation  is 
then  used  to  communicate  nodal  forces  for  the  interface  nodes  to  the  adjacent  subdomains. 
Note  that  the  interfaces  nodes  are  integrated  redundantly  by  all  processors  which  contain 
elements  connected  to  those  nodes.  However,  this  extra  cost  is  compensated  by  the 
reduction  in  communication  costs  brought  about  by  the  fact  that  communication  is  required 
only  once  per  time  step.  The  algorithm  for  a  MIMD  implementation  is  shown  in  Table  1. 

Table  1.  MIMD  implementation  of  exchange  algorithm 


loop  over  subdomains  I  (processors) 

loop  over  elements  e  in  subdomain  I 

gather  nodal  velocity  ve 
compute  rate-of-deformation  by  Eq.(9) 
compute  stress  by  constitutive  equation 
compute  nodal  forces 
assemble  nodal  forces 


end  loop  on  e 
end  loop  on  1 

exchange  forces  on  interface  nodes  between  subdomains 

compute  accelerations  by  Eq.  (5) 

integrate  to  obtain  new  velocities  and  displacements 


Incidentally,  this  algorithm  was  not  very  effective  on  the  CMS,  probably  because  of 
inefficiencies  in  user-directed  message  passing.  However,  on  the  DELTA  this  algorithm 
proved  quite  efficient. 

The  equations  were  integrated  in  time  by  the  explicit  central-difference  method. 
Two  types  of  elements  were  used  in  these  computations; 

1.  Quadrilaterals  with  one-point  quadrature  and  consistent  hourglass  control 
(Flanagan  and  Belytschko(1981)),  with  a  perturbation  stabilization  parameter  e  =  0.01. 

2.  Constant  strain  triangles  arranged  in  the  crossed-diagonal  pattern  which  avoids 
volumetric  locking,  Nagtegaal,  Paries  and  Rice(1974). 

A  viscoplastic  stress-strain  law  in  which 

a  =  -  Tivp)  (11) 

was  used;  is  the  elastic  constitutive  matrix  and  T|'T  the  viscoplastic  strain  rate. 

The  von  Mises  yield  surface  has  been  used  to  determine  the  flow  direction  and  the 
rate  of  viscoplastic  strain  can  be  defined  as  follows 


Ti[jP  =  Ti(o,e)— 
acij 

(12a) 

f(ajc)  =  o  -  k  =  0 

(12b) 

where  effective  stress  a  and  effective  viscoplastic  strain  e  are  given  by 
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(13) 


where  s  u>  the  deviatoric  stress  and 


e  = 


Ti''P:Ti''Pdt 


(14) 


and  the  power  law  used  here  to  determine  the  rate  of  effective  viscoplastic  strain  is 


(15) 


where  a  and  m  are  material  properties  and  g(£)  is  the  hardening  function. 

The  rate  tangent  modulus  scheme,  see  Peirce,  Shih  and  Needleman  (1984),  is  used 
for  the  constitutive  update  for  each  time  step,  so  the  viscoplasdc  strain  is  updated  as 
follows: 

Ai  =  At[(  1  -  0)  Et  +  8£t+At]  ( 1 6) 

where  the  parameter  9  can  vary  from  0  to  1;  the  6  =  0  and  0  =  1  are  the  forward  and 
backward  Euler  method,  respectively.  We  choose  0  =  0.5  here,  which  corresponds  to  a 
semi-implicit  central  difference  update.  To  maintain  accuracy  in  the  later  stages  of 
localization,  the  time  step  should  be  decreased,  because  the  exponential  growth  of  the 
strains  in  the  band  can  lead  to  significant  truncation  errors. 

1.3.  NUMERICAL  SOLUTIONS 

Calculations  were  performed  for  a  square  slab  in  plane  strain  as  shown  in  Fig.  1.  A 
velocity  was  prescribed  at  the  ends  as  shown:  after  a  ramp,  the  constant  value  vq  is  30.0 
m/sec.  The  problem  specification  is  shown  in  Fig.  1;  it  is  similar  to  Needleman  (1989). 
Most  of  the  calculations  were  made  on  a  CM2(X). 

In  all  cases,  an  initial  imperfection  was  included  in  the  form  of  a  reduction  in  the 
yield  strength.  One  of  the  imperfections  is 

aY(x,y)  =  OY[  1  -  ocexp{-  (x  -  xq)^  -  (y  -  yo)^  }/ro]  (17) 
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where  xq,  yo  specify  the  location  of  the  maximum  imperfection  and  Tq  its  radius.  In  the 
problems  reported  here,  xq  =  yo  =  0  except  where  noted  and  a  =  0.20.  The  same 
viscoplastic  material  law  as  reported  by  Needleman  (1989)  was  used;  Young's  modulus  E 
=  211  GPa,  =  460  MPa,  Poisson’s  ratio  n  =  0.3,  the  viscoplastic  exponent  m.  =  0.01. 

In  the  first  example,  two  planes  of  symmetry  were  used  and  only  the  upper  right- 
hand  quadrant  was  modeled.  Unless  noted,  a  value  ro  of  0.1  mm  was  used;  this  will  be 
called  the  large  imperfection.  A  deformed  26  x  20  quadrilateral  element  mesh  is  shown  in 
Fig.  3a  and  exhibits  no  localization;  Needleman  obtained  shear  banding  by  using  a  similar 
refinement  in  a  crossed-triangular  mesh  and  aligning  the  diagonal  of  quadrilateral  elements 
with  the  expected  direction  of  shear  band.  The  results  we  have  obtained  for  the  crossed- 
triangular  mesh,  shown  in  Fig.  3b,  also  show  localization.  Note  that  most  of  the 
deformation  occurs  in  a  single  band,  so  the  width  of  the  shear  band  is  only  slightly  larger 
than  the  element  size.  This  is  a  consequence  of  the  fact  that  the  effective  size  of  the 
imperfection  is  of  the  order  of  the  element  size.  In  these  results,  the  strain  field  is  a  step 
function  in  space  since  the  elements  are  constant  strain  elements  and  one  or  two  elements 
span  the  entire  shear  band. 

The  same  problem  was  then  solved  with  a  128x128  mesh,  i.e.  16K  equally  sized 
elements.  Twofold  symmetry  was  again  used  in  this  problem,  so  only  the  upper  right-hand 
quadrant  is  modeled.  The  evolution  of  the  mesh  for  the  upper  right-hand  quadrant  is 
shown  in  Figs.  4  through  6.  The  figures  show  the  development  of  a  shear  band 
approximately  10  elements  in  width.  Figs.  7  and  8  show  contour  plots  of  the  effective 
viscoplastic  strain,  from  which  it  can  be  seen  that  the  deviation  from  uniform  strain  first 
occurs  over  a  band  of  8  to  10  elements.  The  structure  of  the  strain  field  varies  somewhat 
along  the  length  of  the  band  initially,  but  becomes  quite  uniform  in  the  later  stages,  except 
at  the  two  points  where  the  shear  bands  intersect  (at  the  center  of  the  plate  and  at  the  top 
surface).  With  time,  the  region  of  largest  strain  narrows.  Figure  9  shows  the  deformed 
mesh  for  the  same  problem  with  a  80x80  mesh.  As  can  be  seen  by  comparison  with  Fig. 
4,  the  deformed  meshes  almost  appear  identical. 

Figure  10  shows  the  profiles  of  the  effective  viscoplastic  strain,  given  by  Eq.1.14, 
across  the  shear  band;  this  quantity  gives  a  goood  measure  of  the  severity  of  deformation. 
Results  are  also  shown  for  when:  (1)  each  quadrilateral  element  is  replaced  by  4  constant 
strain  triangular  elements  in  a  cross-diagonal  pattern,  so  the  total  number  of  elements  is 
64K;  (2)  quadrilaterals  using  selective-reduced  integration  are  used.  The  plots  were  made 
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using  the  effective  viscoplastic  strains  in  the  elements  which  lie  along  the  line  Y  =  1 .0  + 
aX,  where  a  is  chosen  so  the  line  is  perpendicular  to  the  direction  of  shear  band.  It  can  be 
seen  that  the  effective  plastic  strain  profiles  for  the  different  elements  are  almost  identical 
for  t  =  2.0iisec  and  for  t  =  3.0|isec.  At  later  times,  a  significant  difference  develops,  and 
the  maximum  effective  strain  differs  markedly  for  the  three  different  elements.  However,  it 
should  be  noted  that  the  growth  of  the  strain  is  exponential  just  as  in  the  one  dimensional 
problem;  see  Belytschko,  Moran  and  Kulkami(1991)  for  a  one-dimensional  solution.  As  a 
consequence,  the  strain-rate  increases  with  the  strain  and  grows  most  rapidly  at  the  point 
where  it  is  maximum.  Therefore,  the  band  appears  to  narrow  in  time  and  the  highest 
gradients  occur  at  the  center,  where  the  strain  profile  becomes  cusp-shaped.  At  4.0p,sec, 
the  cusp-shaped  region  has  narrowed  to  1  or  3  elements.  The  resolution  in  the  cusp-shaped 
region  is  then  insufficient,  and  the  peak  strain  becomes  very  dependent  on  the  type  of 
element.  Thus,  inevitably,  after  a  sufficiently  long  time,  these  meshes  becomes  inadequate 
for  the  viscoplastic  strain-softening  problem. 

The  problem  was  repeated  with  an  80x80  mesh  using  the  same  stabilization 
parameters  (e  =  0.01)  and  a  larger  stabilization  parameter  (e  =  0.05).  The  strain  profiles  are 
shown  in  Fig.  11.  Comparison  of  the  results  with  e  =  0.01  to  those  for  the  128x128  mesh 
(top  of  Fig.  10)  shows  that  the  strain  profiles  for  for  t  =  2.0^sec  and  for  t  =  3.0|isec  agree 
closely.  Similarly,  for  the  two  values  of  the  stabilization  parameter  the  first  two  results  are 
almost  identical.  However,  the  result  at  t  =  4.0^sec  differs  markedly.  This  marked 
difference  is  also  attributed  to  insufficient  resolution  when  the  cusp-shaped  region  only 
spans  a  few  elements. 

The  load  deflection  curves  for  meshes  of  different  size  using  the  one-point 
quadrature  element  with  stabilization  are  compared  in  Fig.  12.  It  can  be  seen  that  the 
results  are  almost  identical  for  the  three  meshes  until  the  latest  stage  of  the  deformation;  the 
horizontal  axis  after  t  =  l.Ojisec  corresponds  to  time  because  the  velocity  is  constant  The 
load  deflection  curves  for  the  one-point  quadrature  element  with  stabilization(URI),  the 
constant  strain  triangles  in  the  cross-diagonal  pattern  and  the  SRI  element  are  compared  in 
Fig.  13.  These  results  are  also  almost  identical  except  for  the  later  stages;  it  is  interesting  to 
note  that  the  results  for  the  triangular  element  and  the  URI  element  are  almost  identical 
throughout  the  simulation. 

The  problem  was  repeated  with  a  128x128  quadrilateral  mesh  with  the  imperfection 
of  Eq.  1.17  with  a  =  0.20,  ro  =  0.0039mm,  i.e.  a  smaller  imperfection.  The  deformed 


9 


mesh  is  shown  in  Fig.  14,  the  effective  viscoplastic  strain  profile  in  Fig.  15.  The  change  in 
the  response  (as  compared  to  Figs.  7  and  8)  is  quite  dramatic:  the  reduction  of  the  size  of 
the  imperfection  reduces  the  width  of  the  shear  band  and  introduces  4  additional 
intersecting  shear  bands  in  the  quadrant.  Apparently,  with  the  smaller  imperfection,  the 
deformation  due  to  the  shear  band  through  the  imperfection  is  not  sufficient  to 
accommodate  the  prescribed  displacement  of  the  top  surface  and  additional  shear  bands  are 
triggered.  This  phenomenon  has  also  been  observed  experimentally,  Anand,  Kim  and 
Shawki(1987). 

A  much  smaller  imperfection,  r©  =  0.{X)195mm,  was  also  used  for  the  same 
problem.  Again,  an  interesting  result  was  found.  Although  the  shear  band  pattern  is  the 
same  as  in  the  previous  results,  the  larger  viscoplastic  strains  now  occur  in  the  shear  bands 
which  do  not  pass  through  the  imperfection;  see  the  contour  plots  of  the  effective 
viscoplastic  strain.  Fig.  16.  This  result  again  shows  the  crucial  role  of  imperfections  in  the 
structure  of  the  shear  bands. 

To  further  assess  the  effects  of  the  size  of  the  imperfection,  the  results  were 
repeated  with  the  imperfections: 

crY(x,y)  =  CTY  { 0.9  +  0.0 1  [(x  -  xq)^  -  (y  -  yo)^  ]/ro }  r<ro  ( 1 8) 

=  ^  r>ro 

with  ro=0.078mm  and  ro=0.156mm;  these  correspond  to  10  and  20  elements,  respectively. 
The  effective  viscoplastic  strain  profiles  at  4.0Msec  are  given  in  Fig.  17.  It  can  be  seen  that 
the  width  of  the  shear  band  is  approximately  twice  as  great  for  ro=  0.156mm.  It  is 
interesting  to  note  that  the  peak  strain  is  equal  for  these  two  cases  at  t  =  4.0|isec;  it  is 
evidently  set  by  the  maximum  magnitude  of  the  imperfection,  which  in  this  case  is 
identical.  The  maximum  magnitude  of  the  imperfection  is  smaller  than  for  the  calculation 
shown  in  Fig.  10,  and  so  the  maximum  strain  is  lower. 

The  location  of  the  imperfection  also  plays  an  important  role  in  the  shear  band 
pattern.  Figure  18  shows  the  shear  band  pattern  in  the  upper  half  of  the  plate  when  an 
imperfection  given  by  Eq.l7  with  a  =  0.2  and  ro=0.10mm  is  placed  at  xo=  -0.25mm, 
yo=0.0mm.  The  shear  band  pattern  now  loses  symmetry  with  respect  to  the  x-plane  in 
order  to  pass  through  the  imperfection.  Similar  patterns  are  observed  when  the 
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imperfection  is  placed  anywhere  on  the  y=0.0  axis.  However,  when  an  imperfection  is  not 
placed  along  the  centerline,  no  shear  bands  are  observed  in  the  computation  for  a  long  time 
(our  computations  were  of  course  limited  in  duration).  For  example,  when  the 
imperfection  of  Eq.(17)  was  placed  at  x=0.1mm,  y=  1.0mm  no  shear  band  pattern  was 
triggered  during  the  first  4.0|isec.  Thus,  the  appearance  and  structure  of  the  shear  band  in 
these  dynamic  problems  depends  on  an  interplay  between  wave  focusing  and  the  position 
of  the  imperfection. 

1.4.  CONCLUSIONS 

The  results  we  have  obtained  show  the  following: 

1.  imperfections  play  an  important  role  in  the  setting  the  external  and  internal 
structure  of  shear  bands  in  viscoplastic  materials  and  can  be  said  to  set  a  length  scale  for 
viscoplastic  materials  in  that  the  dimensions  of  the  imperfection  set  the  initial  width  and 
structure  of  the  shear  band; 

2.  little  difference  is  found  in  the  results  obtained  with  different  elements  or 
different  values  of  the  stabilization  parameter  as  long  as  the  mesh  is  sufficiently  refined  so 
that  more  than  five  elements  occur  within  the  high  gradient  region  of  the  shear  band. 

For  coarse  meshes,  quadrilateral  elements  tend  to  exhibit  less  of  a  tendency  to 
bifurcate  than  cross-diagonal  triangular  meshes;  this  was  also  noted  by  Ortiz,  Leroy  and 
Needleman(1987)  and  Belytschko,  Fish  and  Englemann(1988).  It  has  been  shown  here 
that  for  meshes  in  which  the  elements  are  small  compared  to  the  scale  of  the  imperfection, 
mesh-dependence  for  the  onset  and  early  stages  of  localization  behavior  is  insignificant. 
The  results  obtained  here  show  that  quadrilaterals  with  different  quadrature  schemes  and 
triangular  elements  give  similar  solutions  as  long  as  the  mesh  is  sufficiently  fine  to  capture 
the  high  strain  gradients  in  the  center  of  the  band.  Even  for  the  finest  meshes  used  here, 
this  condition  is  met  only  for  about  3  |i.sec,  which  correspond  to  the  early  stages  of 
localizadon. 

However,  in  the  uncoupled  viscoplatic  localization,  the  strain  profiles  develop  a 
cusp-shaped  profile  which  narrows  with  time.  When  the  cusp  becomes  narrow  enough 
relative  to  element  size,  the  resolution  is  insufficient  and  the  maximum  effective  strain  and 
the  structure  of  the  band  start  to  depend  markedly  on  the  type  of  element,  mesh  refinement 
and  stabilization  procedure.  It  should  be  stressed  that  this  mesh  dependence  is  strictly  a 
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consequence  of  lack  of  resolution.  To  obtain  reasonable  accuracy,  the  portion  of  the  shear 
band  in  which  the  strain  changes  most  rapidly  must  be  spanned  by  4  to  7  constant  strain 
elements  ( just  as  an  elastic  wave  in  a  finite  element  mesh  must  be  spanned  by  4  to  7 
elements  for  reasonable  accuracy).  Since  the  high  gradient  part  of  the  cusp  narrows  with 
time,  any  fixed  mesh  will  eventually  prove  inadequate.  Thus,  there  is  no  evidence  of  any 
pathological  mesh  dependence  in  these  results.  If  the  region  of  highest  strain  gradient  is 
covered  by  only  one  or  two  elements,  the  well  known  manifestations  of  truncation  error 
appear. 


The  running  time  for  the  uniform  16K  mesh  in  these  explicit  calculations  is  about 
1.4|xsec  for  each  element  cycle  on  the  CM200.  The  crossed- triangular  fine  mesh  requires 
about  7  times  as  much  wall  clock  time  as  the  quadrilateral  mesh.  We  would  expect  an 
increase  of  only  a  factor  of  4,  but  evidently  the  increased  memory  requirements  of  crossed- 
triangular  meshes  reduce  the  speed  further. 

The  location  and  size  of  the  imperfection  have  been  shown  to  play  a  critical  role  in 
both  the  internal  shear  band  structure  and  the  arrangement  and  multiplicity  of  shear  bands. 
By  reducing  the  size  of  the  imperfection,  an  interesting  transition  to  multiple  shear  bands 
was  noted.  This  appears  to  arise  from  the  interplay  of  strain  concentrations  caused  by  the 
wave  motion  and  the  effect  of  the  imperfection.  Such  multiple  shear  bands  have  be  noted 
in  experiments  by  Anand  et  al.  (1987).  The  major  role  of  imperfections  in  setting  the 
internal  morphology  of  shear  bands  has  also  been  demonstrated  in  a  two  dimensional, 
dynamic  problem.  In  previous  studies,  because  of  computer  limitations,  two  dimensional 
meshes  usually  used  elements  of  the  same  size  as  the  imperfection,  so  the  importance  of  the 
imperfection  could  not  be  observed. 
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Captions 


Fig.  1.1.  Description  of  model  problem 

Fig.  1.2.  Stress-strain  functions  for  examples;  note  strain-softening  after  effective  plastic 
strain  i  is  about  0.1 

Fig.  1.3a.  Deformed  mesh  at  6.0  psec  with  quadrilateral  elements;  ro  =  0.1mm 

Fig.  1.3b.  Deformed  mesh  at  b.Opsec  with  crossed-triangular  elements  (only  quadrilaterals 
are  shown);  r©  =  0.1mm. 

Fig.  1.4.  Deformed  mesh  for  128x128  quadrilateral  elements  with  a=0.2  at  t  =  4.79psec. 

Fig.  1.5.  Deformed  mesh  for  128x128  quadrilateral  elements  with  a=0.2  at  t  =  b.OOpsec. 

Fig.  1.6.  Deformed  mesh  for  128x128  quadrilateral  elements  with  a=0.2  at  t  =  7.25)isec. 

Fig.  1.7.  Effective  viscoplastic  strain  in  the  problem  of  Fig.  4  through  6  at  t  =  4.0psec; 

128x128  quadrilateral  mesh,  ro  =  0.1mm. 

Fig.  1.8.  Effective  viscoplastic  strain  in  the  problem  of  Fig.  4  through  6  at  t  =  7.25|isec; 
128x128  quadrilateral  mesh,  ro  =  0.1mm. 

Fig.  1.9.  Deformed  mesh  for  80x80  quadrilateral  elements  with  a=0.2  at  t  =  4.79  psec. 

Fig.  1.10.  Effective  viscoplastic  strain  profiles  across  the  shear  band  for  quadrilaterals(top), 
the  128x128  crossed-diagonal  triangular  element  mesh  (middle)  andquadrilaterals  with 
selective-reduced  integration  SRI  (bottom)  at  t=2.0,  3.0,  4.0  psec,  respectively,  r©  = 
0.1mm;  all  meshes  are  128x128. 

Fig.  1.1 1.  Effective  viscoplastic  strain  profiles  across  the  shear  band  for  quadrilateral  80x80 
with  stabilization  parameters  of  0.01  (top)  and  0.05  (middle). 

Fig.  1.1 2.  Average  stress  versus  end  displacement  curves  for  three  types  of  elements  in  a 
128x128  mesh:  constant  strain  triangles,  quadrilaterals  with  selective  reduced  integration 
(SRI)  and  quadrilaterals  with  one-point  underintegration  and  stabilization(URl). 
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Fig.  1.13.  Average  stress  versus  end  displacement  curves  for  three  meshes  of  one-point 
quadrature  elements;  6,400  elements,  10.000  elements  and  16,384  elements. 

Fig.  1.14.  Deformed  mesh  for  small  imperfection,  ro=0.0039mm,  at  t=7.25  ^sec,  showing 
of  multiple  shear  bands. 

Fig.  1.1 5.  Effective  viscoplastic  strain  for  small  imperfection  at  t=7.25  ^sec; 
ro=0.0039mm. 

Fig.  1.16.  Effective  viscoplastic  strain  across  the  shear  band  for  small  imperfection  at 
t=7.25  ^sec  (ro=0.00 195mm)  showing  the  development  of  multiple  shear  bands  with 
markedly  different  strain  distributions  than  in  Fig.  1.15. 

Fig.  1.17.  Effective  viscoplastic  strain  profiles  across  the  band  for  the  imperfection  Eq.(18) 
at  t=4.0|isec. 

Fig  1.18.  Deformed  mesh  for  upper  half  of  specimen  at  t=5.0^isec  for  imperfection  at  x=- 
0.25mm,  y=0.0;  ro=0.1mm. 
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Fig.  1.3a.  Deformed  mesh  at  6.0  jisec  with  quadrilateral  elements;  ro  =  0.1mm 


Fig.  1.3b.  Deformed  mesh  at  6.0iJ.sec  with  crossed-triangular  elements  (orUy 
quadrilaterals  are  shown);  ro  =  0.1mm. 


Fig.  1.6.  Deformed  mesh  for  125x123  quadrilateral  elements  with  a=0.2  at  t 


Fig.  1.7.  Effective  viscoplastic  strain  in  the  problem  of  Fig.  4  through  6  at  t 
4.0psec;  123x128  quadrilateral  mesh.  To  =  O.lmm. 


Fig.  1.9.  Deformed  mesh  for  SvJxiO  cuadri'ateral  elements  with  a=0.2  at  t  = 
4.79  |isec. 


Fig.  1.10.  Effective  viscoplastic  strain  profiles  across  the  shear  band  for 
quadrilaterals(top),  the  128x128  crossed-diagonal  triangular  element  mesh 
(middle)  andquadrilaterals  with  selective-reduced  integration  SRI  (bottom)  at 
t=2.0, 3.0, 4.0  jisec,  respectively,  ro  =  0.1mm;  all  meshes  are  128x128. 
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Fig.  1.12.  Average  stress  versus  end  displacement  curves  for  three  types  of 
elements  in  a  128x128  mesh:  constant  strain  triangles,  quadrilaterals  with 
selective  reduced  integration  (SRI)  and  quadrilaterals  with  one-point 
underintegration  and  stabili2ation(URI). 
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Fig.  1.13.  Average  stress  versus  end  displacement  curves  for  three  meshes  of 
one-point  quadrature  elements;  6,400  elements,  10,000  elements  and  16,384 
elements. 


Fig.  1.14.  Deformed  mesh  for  sir^all  imperfection,  ro=0.0039mm,  at  t=7.25  )j.sec, 
showing  of  multiple  shear  bands. 


Fig.  1.15.  Effective  viscoplastic  strain  for  small  imperfection  at  t-7.25  jisec 
ro=0. 0039mm. 


Fig.  1.16.  Effective  viscoplastic  strain  across  the  shear  band  for  small 

imperfection  at  t=7.25  psec  (ro=0. 00195mm)  showing  the  development  of 
multiple  shear  bands  with  markedlv  different  strain  distributions  than  in 
Fig.  15. 
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Fig.  1  17.  Effective  viscopiaslic  strain  profiles  across  the  band  for  the 
imperfection  Eq.(18)  at  t=4.0usec. 


CHAPTER  2 


H-ADAPTIVE  FINITE  ELEMENT  METHODS  FOR 

DYNAMIC  PROBLEMS,  WITH  EMPHASIS  ON  LOCALIZATION 

2.1.  INTRODUCTION 

The  aim  of  this  woric  is  to  develop  methods  for  h-adaptive  solutions  of  nonlinear 
transient  problems  with  localization  behavior.  Localization  occurs  in  rate  independent 
materials  when  the  momentum  equations  lose  hyperbolicity;  a  regularization  procedure 
based  on  damping  known  as  a  viscoplastic  model  will  be  used  here  which  prevents  loss 
of  hyperbolicity.  Even  with  regularization,  the  effect  of  a  negative  value  of  modulus  is 
that  severe  deformations  occur  over  narrow  regions  in  which  the  strains  and  displacement 
gradients  are  very  large;  material  response  with  negative  tangent  modulus  are  called 
strain-softening  in  engineering.  Physically  localization  results  in  bands  of  high  strain 
which  are  known  as  shear  bands  in  the  deformation  of  continua  and  as  hingelines  in  the 
deformations  of  shells. 

In  developing  these  methods  we  have  focused  our  efforts  on  explicit  time 
integration  methodologies  and  finite  elements.  Our  focus  on  explicit  time  integration  and 
finite  elements  is  motivated  by  the  fact  that  large  scale  engineering  problems  are  almost 
exclusively  solved  by  such  methods;  DYNA3D  by  Hallquist  and  Benson  (1986)  and  its 
descendants  today  are  the  dominant  tools  for  nonlinear  transient  problems  of  solid 
mechanics  problems. 

There  is  already  a  rich  literature  in  adaptive  methods  for  transient  problems,  so 
only  a  few  works  will  be  noted.  Adaptive  methods  for  hyperbolic  problems  were 
developed  by  Berger  and  Oliger  (1984).  They  solved  one  and  two  dimensional  transient 
problems  with  finite  difference  methods  and  used  an  error  criterion  based  on  Richardson 
extrapolation.  Adjerid  and  Flaherty  (1986)  applied  h-adaptive  finite  difference  methods  to 
parabolic  problems.  Drew  and  Flaherty  (1984)  applied  adaptive  space-time  elements  to 
one-dimensional  shear  banding. 

In  the  finite  element  literature  of  solid  mechanics,  the  application  of  h-adaptive 
methods  to  transient  problems  has  been  sparse,  and  we  know  of  no  applications  of  h- 
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adaptive  finite  element  methods  to  transient  problems  with  unstable  material  behavior,  i.e. 
strain  softening. 

There  have  however,  been  a  number  of  interesting  related  works.  Ortiz  and 
Quigley  (1991)  have  used  an  r-adaptive  method  with  an  error  criterion  based  on 
interpolation,  Diaz  et  al  (1983),  to  solve  some  static  localization  problems.  Zeng  and 
Wiberg  (1992)  developed  an  h-adaptive  method  based  on  constant  strain  triangular 
elements  using  a  modification  of  the  Zienkiewicz-Zhu  (1987)  or  ZZ  error  criterion  for 
linear  wave  problems.  It  has  been  proven  by  Rank  and  Zienkiewicz  (1987)  that  for  linear 
static  problems  this  error  criterion  is  equivalent  to  the  Babuska-Rheinboldt  (1978) 
criterion  of  the  residual  type. 

One  of  the  aims  of  this  paper  is  to  study  various  error  criteria  for  this  class  of 
problems.  The  following  error  criteria  were  selected:  The  Zienkiewicz-Zhu  (1987)  L2- 
projection  on  stress,  the  Babuska-Rheinboldt  (1978)  residual  criterion,  and  the 
interpolation  criterion  Diaz  et  al  (1983).  Though  the  first  two  of  these  were  proposed  for 
linear,  elliptic  problems,  they  appear  to  be  reasonably  effective  in  linear  hyperbolic,  one¬ 
dimensional  problems.  However,  in  localization  problems  they  do  not  fare  as  well.  Our 
studies  show  that  the  ZZ  criterion,  when  formulated  in  terms  of  stresses  is  quite 
ineffective  for  localization  problems.  In  fact,  the  stress  field  in  and  adjacent  to  a 
localization  band  which  is  triggered  by  strain-softening  is  nearly  constant,  so  that  the  ZZ 
criterion  never  initiates  refinement  in  the  area  of  localization.  Residual  type  criteria  are 
similarly  ineffective  in  localization. 

In  this  paper  we  also  propose  an  error  criterion  based  on  a  L2  projection  of  the 
strains.  It  is  shown  that  this  technique  effectively  captures  the  region  where  refinement, 
or  what  we  call  fission  of  elements,  is  needed.  Results  are  obtained  for  several 
localization  problems  and  compared  to  results  obtained  with  very  fine  meshes.  These 
results  show  that  the  present  procedure  is  capable  of  achieving  solutions  of  comparable 
accuracy  with  far  fewer  degrees  of  freedom. 

The  L2  projection  is  more  easily  applied  in  a  general  finite  element  setting  than 
interpolation  criteria;  the  latter  requires  a  knowledge  of  adjacent  elements.  It  also  appears 
to  us  simpler  to  apply  than  the  local  evolution  equation  approaches  pioneered  by  Rank 
and  further  developed  by  Oden  and  coworkers  (1989).  In  the  latter,  it  is  necessary  to 
solve  the  governing  equations  with  a  bubble  mode  added  to  the  element  and  an  estimate  of 
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the  tractions  at  the  element  interface.  Although  these  methods  have  proven  very 
successful  in  two  dimensional  meshes  (and  are  probably  equally  effective  in  three 
dimensional  meshes),  their  application  to  shell  structures  looks  formidable.  It  is  not  clear 
to  us,  for  example,  how  these  methods  would  be  applied  to  a  T-section  or  a  box  beam 
modeled  by  shell  elements,  whereas  the  L2  strain  projection  is  readily  applied  as 
described  in  Belytschko  and  Yeh  (1992).  Although  general  structures  are  not  addressed  in 
this  paper,  the  objective  is  to  develop  error  criteria  suitable  for  transient  analysis  of 
nonlinear  solids  and  continua. 

An  outline  of  the  paper  is  as  follows.  In  Section  2,  the  governing  equations  are 
reviewed.  In  Section  3,  the  error  criteria  used  in  this  study  are  presented  along  with  the 
new  error  criterion  based  on  an  L2  strain  projection.  The  adaptive  technique  is  covered  in 
Section  4.  Numerical  results  for  one  dimensional  and  two  dimensional  problems  are 
presented  in  Section  5,  followed  by  a  discussion  and  conclusions  in  Section  6. 

2.2.  GOVERNING  EQUATIONS 

The  problems  were  solved  by  a  finite  element  method  with  a  Lagrangian  mesh. 
The  rate-of-deformation  tensor 


2  3xj  d\ 


or  T|  =  D  V 


(1) 


is  used  as  a  measure  of  deformation;  here  Vj  are  the  velocities,  Xj  the  Eulerian  (spatial) 

coordinates.  The  material  is  characterized  as  an  elastic-viscoplastic  von  Mises  solid,  the 
constitutive  equation  is  expressed  in  terms  of  the  Jaumann  rate  of  the  Cauchy  stress 


a  =  C(Ti  -  qP) 


(2) 


where  C  is  the  initial  elastic  constitutive  matrix.  qP  is  the  viscoplastic  strain  rate  given  by 


qP  = 


(3) 
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with  the  stress  de viator  o'  =  o  -  ^  a:I,  where  I  is  the  identity  tensor.  The  radius  of  the 
3  1 

flow  potential  is  O  =  (2  o':  <s')  .  Here,  the  rate  of  the  effective  plastic  strain,  t  is 

speciHed  by  the  power  law  relation 


•  •  r-/ 

e  =  Eq  [a/g(e)] 


(4) 


where  g(E)  is  the  hardness  function,  m  is  the  strain-rate  hardening  exponent  and  Eq  is  a 
reference  strain  rate. 

The  momentum  equation  are  Oy  j+  bi  =  pvi  (5a) 

in  rectilinear  coordinates,  it  can  be  written  as 


D  ^o  +  b  =  pv 


(5b) 


where  b  is  the  vector  of  body  forces  and  r  is  the  mass  density. 

The  finite  element  equations  are  obtained  by  approximating  the  velocities  in  each 
element  by; 


v(x,t)  =  N(x)d(t) 


(6) 


The  resulting  ordinary  differential  equations  are: 
momentum  equation 

Ma(t)  =  f(t) 


^(0  =  *ext(0  -  fint(0 


^xt  —  ^  Lj^xt(0 


*in 


t(t)=( 

JQe 


B^(x,t)o(x,t)  dfl 


•  =  §^=d(t) 


fint=SLjf?o,(t) 


(7) 

(8) 

(9) 

(10) 
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measure  of  d^ormation 

Ti(x,t)  =  B(x)de(t)  B=DN  (11) 

de(t)  =  Led(t)  (12) 

In  the  above: 

t  =  time; 

M  =  mass  matrix; 
d  =  nodal  displacement  matrix; 

^t).  (ext(0.  ^mt(0  =  resultant,  external  and  internal  nodal  forces; 

B(x)  =  semidiscrete  form  of  the  symmetric  part  of  the  gradient  operator; 

T]  =  rate-of-deformation  (velocity  strain)  tensor, 

Le  =  Boolean  connectivity  matrix  of  element  e. 

2.3.  ERROR  MEASURES 

The  error  is  thought  to  be  due  to  spatial  discretization  only,  no  attempt  is  done  to 
quantify  time  integration  error.  However,  the  time  integration  scheme  used  here  is  an 
explicit  central  difference  method  which  restricts  the  time  step  to  a  small  critical  value. 
The  error  is  calculated  using  numerical  integration  with  three  point  quadrature  for  the  one 
dimensional  problem  and  2x2  quadrature  for  the  two  dimensional  problems.  We  consider 
four  types  of  element  error  measures: 

(1)  Zienkiewicz-Zhu  (ZZ) 

\ 

e“=  {  J(?y-ol»)’^(?y-oh)dfl}/Qi  ^  (13) 

L  «i 
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This  error  criterion  was  first  proposed  for  linear  elastic  problems  (in  the  original  ZZ 
criterion,  the  elastic  coefficient  matrix  is  included  to  give  an  error  in  energy,  but  it  is 
omitted  here  to  simplify  application  to  nonlinear  problems).  Here  Qi  is  the  volume  of 
element  i,  are  the  finite  element  solution  for  the  stresses,  and  is  the  L2  projection  of 
the  finite  element  stress  which  is  given  by  (see  Oden  and  Brauchli  (1971),  Zienkiewicz 
and  Zhu  (1987)): 

s(x)=  N  s 

where  s  is  the  array  of  the  nodal  value  of  the  projected  stress  "5,  which  is  a  C®  function 
whereas  is  a  C*^  function  in  these  elements. 

s=lVl-ljNTohdO  (15) 

O 


jNTNdO  (16) 

Q 

Rl  is  similar  to  the  mass  matrix  and  a  lumped,  diagonal  form  of  the  mass  matrix  obtained 
by  Lobatto  quadrature  was  used. 


(2)  Babuska-Rheinboldt  (BR) 


This  is  an  error  criterion  developed  for  linear  elliptic  (static)  problems.  In  applying  their 
concept  to  transient  problems,  we  consider  the  residual  in  the  equations  of  motion. 


{Cl  Jr2dO+c2  Jj2dr}/ni 

«i  n 


(17) 


where  in  one  dimensional  problems 


‘^l~24k  ‘^2  =  24k  ^ 


1-v 


3i|- 


r  =  a,x  +  b  -  pv 


J  =  <  r  > 

where  hi  is  the  length  of  element  i,  n  is  Poisson's  ratio  and  <>  is  the  jump  operator,  i.e. 
J  is  the  inter-element  traction  jumps. 


(3)  Interpolation  (IN) 

This  criterion  was  developed  by  Diaz  et  al  (1983)  for  static  problems.  The  form  we  use 
here  for  dynamic  problems  is  identical  and  in  one  dimension  is  given  by: 


m 
e . 

1 


i  f(hiU*xx)^dx}/hi 

hi 


1_ 

|2 


(18) 


where  u,^jj  is  obtained  by  using  a  quadratic  Lagrange  interpolant  for  constant  strain 

elements.  In  general,  a  Lagrange  interpolant  one  order  higher  than  that  used  for  the  finite 
element  solution  is  used.  In  one  dimension 

2ui  2u2  2u3 

“,xx  =  (xi-X2)(xi-X3)  (x2-xi)(X2-X3)  ^  (x3-xi)(x3-X2) 

xj,  X2,  and  X3  are  the  coordinates  of  midpoints  of  elements  i-1,  i,  i+1  for  an  interior 
element  i.  For  an  element  i  with  a  boundary  on  the  right-hand  side,  the  midpoints  of  i-2, 
i-1,  i  are  used.  For  an  element  i  with  a  boundary  on  the  left-hand  side  the  midpoints  of  i, 
i+1,  i+2  are  used.  The  displacements  at  xj,  X2,  and  X3  are  denoted  uj,  U2,  U3. 


(4)  Strain>Projection  (SP) 
The  criterion  proposed  here  is 


j(f-eh)V-eh)dn}/Qi 


a 


(20) 


where  F  is  evaluated  by  an  L2  projection  which  is  identical  to  the  L2  projection  ofS.  This 
measure  represents  the  R.M.S  strain  error  in  an  element  and  is  convenient  to  use  because 
it  is  nondimensional.  The  strain  projection  cannot  be  used  directly  at  material  interfaces 
because  strains  are  generally  discontinuous  at  such  interfaces.  The  L2  projection  must 
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then  be  computed  separately  in  each  material  subdomain;  this  procedure  has  been 
described  by  Belytschko  and  Yeh  (1992). 

2.4.  ADAPTIVE  PROCEDURE 

We  use  an  h-adaptive  technique,  we  shall  refer  to  the  refinement  of  an  element  as 
■fission’  and  to  the  unrefinement  of  a  group  (family)  of  four  elements  (siblings)  to  their 
parent  element  as  'fusion',  see  Belytschko  et  al  (1989).  These  processes  are  illustrated  in 
Fig.  la.  When  an  element  is  fissioned  next  to  an  unfissioned  element,  slave  nodes  are 
created.  The  motion  of  a  slave  node  A  is  governed  by  the  constraint  of  compatibility,  i.e. 

Va  =  t{^‘}  (21) 

where  T  is  a  linear  operator  which  enforces  compatibility  and  V  i  and  V2  are  the 
velocities  of  the  master  nodes.  When  node  A  is  midway  between  nodes  1  and  2,  T  =  [  ^  I 

,  ^  I].  The  equations  of  motion  are  not  evaluated  at  slave  node  but  instead  the  nodal 
forces  at  slave  nodes  are  added  to  forces  at  the  corresponding  master  nodes,  i.e. 

(2.22) 

where  fA  are  the  nodal  forces  at  the  slave  node  and  f*  are  the  nodal  forces  at  nodes  1  and 
2  prior  to  the  consideration  of  fA-  This  is  a  standard  technique  for  treating  constraints  in 
explicit  methods. 

The  adaptive  process  is  driven  by  element  error  measures.  An  element  is  a 
candidate  for  fission  when  the  error  is  more  than  a  specified  fission-limit  and  a 
contiguous  group  of  elemr  .its  originating  from  a  single  parent  element  is  a  candidate  for 
fusion  when  the  error  in  all  its  siblings  is  lower  than  a  specified  fusion-limit  e^“.  A 
candidate  is  only  fissioned  if  it  complies  with  the  1-irregular  rule  (Devloo  et  al  (1987)) 
and  the  hmin  requirement,  where  hmin  is  a  prescribed  lower  limit  on  the  size  of  an 
element.  A  candidate  is  fused  if  it  complies  with  the  1 -irregular  rule  only. 

An  important  ingredient  of  an  h-adaptive  program  is  the  data  structure.  In 
designing  a  data  structure,  compromises  must  be  made  to  conserve  storage  and  avoid 
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excessive  computation  time.  In  order  to  simplify  the  data  structure  we  will  enforce  the  1- 
irregular  rule,  i.e.  neighboring  elements  can  be  at  most  one  generation  apart.  The  nodal 
and  element  data  for  the  adaptive  process  is  handled  by  three  arrays: 

(1)  NSBLNG(4,NF):  the  four  siblings  of  a  family  (Fig.  la); 

(2)  NABOR(8,NE);  the  eight  neighbors  of  an  element  (Fig.  lb); 

(3)  MGEN(NE);  the  number  of  descendant  generations  for  an  element. 

where  NE  and  NF  are  the  number  of  elements  and  families  in  the  mesh.  NF  is  initially  set 
to  zero  and  increases  (decreases)  by  one  for  each  fissioned  element  (fused  family). 
NSBLNG  is  used  in  the  fusion  process  in  order  to  identify  elements  of  a  family.  NABOR 
is  used  in  both  fission  and  fusion  to  determine  the  appropriate  numbering,  connectivity 
and  status  of  nodes;  it  is  also  used  to  identify  nodes  lying  on  a  boundary.  MGEN  is 
required  in  order  to  comply  with  the  1 -irregular  rule. 

In  order  to  clarify  the  use  of  the  above  arrays,  let  us  consider  mesh  'a'  shown  in 
Fig.  Ic  with  element  5  already  fissioned.  Also  consider  mesh  'b'  obtained  by  fissioning 
element  2,  The  NABOR  array  for  element  2  is  given  in  Table  1  for  both  meshes.  The  first 
thing  to  notice  is  the  zero  entries  for  L7  and  L8,  indicating  a  boundary  to  the  left-hand 
side  of  element  2;  thus  a  new  node  needs  to  be  created  (node  20),  which  is  a  master  node. 
Next,  L1=L2=1  indicating  that  element  1  extends  along  the  entire  side  of  element  2  and 
thus  a  new  node  (node  18)  needs  to  be  created,  which  is  a  slave  node.  Now,  L3=5  and 
L4=9  thus  a  node  already  exists  (node  16)  which  was  a  slave  node  in  mesh  'a'  and  will 
become  a  master  node  upon  fission  of  element  2.  The  NSBLNG  array  (Table  2)  gives  the 
four  siblings  for  the  only  family  in  mesh  'a',  so  elements  5,7,8  and  9  are  the  only 

n 

members  of  the  mesh  whose  error  is  checked  against  e^^. 

We  will  implement  the  repeat  algorithm  employed  by  Belytschko  and  Yeh  (1992). 
The  basic  idea  is  as  follows.  Start  at  time  t=0  with  an  initial  coarse  mesh  'a'  (Fig.  2)  and 
solve  the  problem  using  a  time  increment  Dt  for  an  elapsed  time  ti,  during  which  the 
mesh  is  modified  at  times  at^,  a=l/n,  n>l;  and  ends  up  as  mesh  'b'.  Now,  we  go  back  to 

t=0  and  use  mesh  'b'  to  solve  the  problem  and  record  the  results.  So  the  first  run  over  the 
time  interval  serves  the  purpose  of  finding  a  suitable  mesh  while  the  second  run  is  used  to 
actually  solve  the  problem  at  hand.  This  procedure  is  repeated  over  every  lime  interval  ti 
as  illustrated  in  Fig.  2  thus  using  mesh  'c*  to  go  from  t=ti  to  t=2ti;  and  so  on  until  the 
end  time  t2  is  reached. 
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Although  the  above  description  was  limited  to  a  four-node  quadrilateral,  it  is  also 
applied  to  the  two  node  element  for  the  one  dimensional  problems  solved  here.  Fig.  la 
shows  the  fission-fusion  process  for  one  element,  in  this  case  there  is  no  need  for  slave 
nodes,  i.e.  all  nodes  are  masters.  Fig.  lb  shows  LI  and  L2  the  two  neighbors  of  element 
JE.  The  same  three  arrays:  NSBLNG,  NABOR,  MGEN  are  used,  however  the  size  of 
NSBLNG  reduces  to  2xNF  and  that  of  NABOR  to  2xNE. 

2.5.  NUMERICAL  EXAMPLES 

Stability  considerations 

The  explicit  integration  schemes  for  the  semi-discrete  equations  of  motion  (Eqn. 
7)  and  the  plastic  strain  rate  (Eqn.  4)  are  conditionally  stable.  In  addition  since  h- 
adaptivity  introduces  much  smaller  elements,  subcycling  in  which  different  time  steps  are 
used  in  different  elements  was  used  to  achieve  computational  efficiency,  see  Belytschko 
and  Liu  (1982)  and  Belytschko  and  Lu  (1992). 

One  dimensional  bar 

Consider  a  bar  of  length  L=100  in,  cross-sectional  area  A=1.0  in^,  elastic 
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modulus  E=:3xl0^  lbf(in2,  density  r=0.00074  Ib/in^.  The  element  mass  matrix  is  M  =  ^  ( 

Mdiag  +  ivicons)  where  M^iag  the  diagonal,  lumped  mass  and  M^ons  the  consistent 
mass.  It  is  well  known  that  this  combination  provides  the  best  spectral  fidelity  for  the 
mesh. 


The  material  is  rate-dependent  with  an  effective  plastic  strain  rate  represented  by 
the  power  law  relation  given  by  Eqn.  4.  The  evolution  equation  (Needleman  1988)  for 

g(E)  takes  the  form  g  =  H  lei,  where  H  =  Hj  for  e  <  e^jj,  H  =  H2  for  e  >  and  g(0)  = 
Oq.  At  each  time  increment  the  stress  state  is  updated  using  the  forward  Euler  scheme. 
Fig.  3  shows  computed  stress-strain  curves  (A,  B,  and  C)  at  various  imposed  strain  rates 
for  Cq=  0.()02E,  m=0.02,  eQ=0.002,  Hj=25aQ,  H2=2.5aQ  and  ejjj=0.2aQ/E.  As 

apparent  from  these  curves  the  material  hardens  with  increasing  strain,  but  because  of 
strain-rate  dependence  softening  might  still  occur. 

The  bar  is  subjected  to  a  constant  velocity  vo=300  in/s  at  both  ends  resulting  in 
two  step  waves  of  constant  strain  traveling  toward  each  other.  The  magnitude  of  the 
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velocity  is  chosen  so  that  the  material  remains  elastic  until  the  waves  meet  at  the  middle  of 
the  bar,  doubling  the  strain  value  at  that  point  and  driving  the  material  well  into  the  plastic 
range.  At  first  thought,  no  localization  is  expected  since  the  material  is  hardening  under 
constant  strain  rates,  but  a  closer  look  at  the  actual  stress-strain  response  (Curve  D,  Fig. 
3)  for  the  midpoint  under  this  particular  loading  shows  a  softening  behavior  which  is  due 
to  the  step  strain  wave.  Solving  this  example  using  a  fixed  mesh  with  80  elements, 
At=2.SxlO*^,  and  t2/At=18(X)  results  in  a  localized  strain  profile  B  in  Hg.  4a  and  a  stress 
profile  as  shown  in  Fig.  4b.  We  will  use  these  results  as  a  benchmark  for  the  adaptive 
results.  These  profiles  are  plotted  for  t=0.9to,  where  to  is  the  time  required  for  an  elastic 
wave  to  traverse  the  entire  length  of  the  bar,  to=L/VE7r . 

It  is  important,  however,  to  keep  in  mind  that  in  problems  involving  localization 
due  to  material  softening,  the  late  stages  of  the  solution  depend  on  the  element  size 
employed  in  the  analysis.  A  simple  remedy  for  this  is  to  restrict  the  minimum  size  of 
elements  used  to  a  minimum  allowable  size;  we  used  hniin=L/80  for  all  runs.  This  is 
appropriate  here,  since  our  goal  is  to  illustrate  adaptivity  and  not  to  propose  a  constitutive 
law  that  overcomes  this  dependence  on  mesh  size  which  occurs  in  the  late  stages  of 
localization  even  with  viscoplasticity. 

In  the  results  for  the  fixed  mesh  all  four  estimates  indicate  high  error  at  the  wave 
fronts  before  the  two  waves  meet  (t<0.5to)-  However,  after  the  waves  meet.  Fig.  5,  only 
the  SP  and  IN  indicate  high  errors  at  wave  fronts  and  localization  regions,  while  the  ZZ 
and  BR  indicate  high  errors  at  wave  fronts  only,  thus  missing  the  localization.  Although 
the  above  error  profiles  are  shown  for  t=0.9to,  they  are  representative  of  profiles  for 
t>0.5to.  The  SP  and  IN  estimates  are  normalized  with  respect  to  Eq  =  while  the  ZZ 

and  BR  are  normalized  with  respect  to  Sq.  The  normalized  SP  estimates  the  error  at  wave 
fronts  to  be  around  0.05  (Fig.  5a)  which  is  the  same  as  the  estimate  given  by  the 
normalized  ZZ  (Fig.  5c),  this  is  expected  because  in  this  region  the  behavior  is  elastic  and 
the  two  estimates  are  identical.  The  normalized  ZZ  and  BR  profiles  (Fig.  5c,d)  are  very 
close,  which  is  expected  in  linear  static  problems  and  seems  to  carry  over  to  dynamics  in 
this  speciHc  example. 

This  one  dimensional  example  was  also  solved  using  an  adaptive  mesh  guided  by 
ZZ  (e^^=0.02a5  and  e^^=0.5e^^)  and  SP  (e^^=0.02eo  with  At=2.5xl0-7, 

ti/Dt=100,  t2/At=1800,  a=0.2,  and  hniin=  L/80.  The  evolution  of  the  adaptive  meshes 
with  time  are  shown  in  Fig.  6.  The  meshes  are  similar  until  localization  occurs  (t=0.5to), 

39 


the  ZZ  driven  mesh  (Fig.  6a)  does  not  reHne  the  mesh  to  accommodate  the  localization 
while  the  SP  does  (Fig.  6b).  Thus  the  strain  profiles  after  localization  are  not  the  same, 
with  the  SP  driven  one  being  closer  to  the  fixed  mesh  (Fig.  4a).  Results  obtained  by 
using  IN  are  similar  to  that  of  SP,  while  results  obtained  by  using  BR  are  similar  to  that 
ofZZ. 


Compression  of  a  rectangular  block 

Consider  the  problem  of  plane  strain  compression  analyzed  by  Needleman  (1989) 
for  dynamic  shear  band  development  The  main  difference  between  Needleman's  analysis 
and  ours  lies  in  the  type  of  discretization;  we  use  four-node  quadrilaterals  with  one  point 
integration  and  hour-glass  control  (Flanagan  and  Belytschko  (1981)),  Needleman  used 
crossed  triangular  quadrilaterals.  We  use  one  point  integration  and  perturbation  hourglass 
control  (with  a  viscosity  factor  of  0.01)  in  order  to  be  able  to  represent  incompressible 
deformation  without  mesh  locking  and  for  speed  of  calculation. 

The  uniaxial  response  for  the  effective  plastic  strain  rate  is  represented  by  a  power 
law  rate  relation  (Eqn.  4)  with  a  hardness  function 


g(e)  =  <^0 


(1  +e/eQ ) 


N 


1  +  (g/ep-^ 

An  initial  imperfection  is  included  in  the  form  of  a  reduction  in  the  yield  strength 


(23) 


(24) 


The  material  properties  used  are  as  follows:  Young’s  Modulus  E=211GPa,  Poisson's 
ratio  v=0.3,  mass  density  p  =  7800  kg/m^,  y=0.2,  ro  =0.1mm,  <jQ=460MPa, 

£^=0.002 18,  N=0.1,  m=0.01  and  £^=0.002  s*!. 

We  analyze  one  fourth  of  the  block  (hoxho)  shown  in  Fig.  7  by  making  use  of 
symmetry.  The  prescribed  velocity  v(t)  is  equal  to  vot/t©  for  t<to  and  vq  for  t>to  with 
to=0.1xl0‘^.  The  end  displacement  and  the  corresponding  average  stress  are  given  by: 

t  .  ^0 

u  =  J  v(t)  dt  and  Cavg  ~  h  J 


(25) 


where  ty  is  the  traction  in  the  y-direction  and  h  is  the  current  width  at  y=ho.  At  each  time 
increment,  the  stress  state  is  updated  using  the  rate  tangent  method  of  Peirce  et  al.  (1984) 
with  the  integration  parameter  q=1.0. 

Case  1 

As  a  first  step,  we  used  a  fixed  mesh  of  20x20  elements  with  two  forms  of  the 
hardness  function  g(£);  (a)  hardening  with  ei=10^EQ  and  (b)  softening  with  ei=102eQ  . 

Fig.  8  shows  average  stress  versus  end  displacement,  with  an  imposed  velocity  vo=10 
m/s,  At=2.0xl0‘9,  t2/At=6000  for  the  two  materials  (Curves  A,C).  The  response  for 
both  materials  is  initially  identical;  however,  for  the  softening  material  there  is  an  abrupt 
stress  drop  with  increased  loading  that  is  associated  with  shear  band  formation.  For  the 
softening  solid  a  region  of  high  effective  strain,  £,  starts  developing  at  -u/ho=0.06  and 
acquires  the  shape  shown  in  Fig.  9a  for  the  last  calculated  point,  -u/ho=0.12.  A  definite 
band  with  high  strains  forms  at  45°  to  the  compression  axis.  The  strains  inside  the  band 
increase  significantly  with  continued  loading,  while  regions  outside  this  band  experience 
little  change  in  strain.  This  is  not  the  case  for  the  hardening  solid,  where  strains  increase 
uniformly  throughout  the  solid  with  increased  loading  (Fig.  9b).  However,  only  a  small 
region  of  higher  strains  develops  near  the  inhomogeneity. 

Fig. .  10  shows  the  surfaces  of  the  SP  error  corresponding  to  the  surfaces  of 
effective  strains  shown  in  Fig.  9.  For  the  softening  solid,  the  error  is  high  in  the  entire 
shear  band  region  (Fig.  10a).  For  the  hardening  solid,  the  SP  error  is  large  only  at  the 
inhomogeneity  location  (Fig.  10b).  These  results  are  desirable  and  will  lead  to  a  useful 
mesh:  refinement  occurs  near  the  inhomogeneity  for  hardening,  and  within  the  whole 
shear  band  for  softening. 

Starting  with  an  initial  mesh  of  20x20  elements,  and  using  the  SP  error  measure 
to  guide  adaptivity  with  e^‘=  0.01,  6^**=  0,  Dt=1.0xl0’9,  ti/At=800,  t2/At=8000,  and 
ot=0.5  the  softening  solid  was  solved.  The  response  curve  D  is  shown  in  Fig.  8.  Labels  1 
to  6  are  used  to  show  the  evolution  of  the  mesh  with  increased  loading.  Each  segment  of 
the  curve  between  labels  represents  a  different  mesh.  For  the  first  segment  (0-1)  no 
refinement  is  needed,  for  segment  1-2  only  few  element  at  the  lower  left  comer  of  the 
mesh  are  fissioned.  For  segment  2-3  the  refinement  is  minimal  and  is  limited  to  a  region 
near  the  imperfection  (Fig.  1  la).  However,  for  segments  3-4, 4-5,  and  5-6  which  involve 


localization,  the  entire  shear  band  region  is  refined  to  give  the  meshes  in  Fig.  1  lb,  c,  d 
respectively.  In  this  analysis  we  use  hniin=ho/80,  which  allows  for  only  two  levels  of 
fission. 


Curves  B,  C  and  D  (Fig.  8)  are  for  the  same  problem,  the  only  difference  being  in 
the  size  of  the  smallest  elements  in  the  shear  band.  Although  the  curves  are  identical  up  to 
about  -u/ho=0.05,  the  response  becomes  noticeably  different  as  the  loading  progresses. 
In  order  to  compare  the  responses  obtained  for  -u/ho>0.05  we  consider:  (1)  the  onset  of 
localization  to  be  defmed  as  the  point  on  the  curve  after  which  the  Hrst  abrupt  drop  in 
stress  occurs,  and  (2)  the  magnitude  of  this  drop  to  be  a  gauge  for  the  severity  of  the 
localization.  Curve  D  predicts  an  earlier  onset  of  localization  with  high  severity  followed 
by  high  oscillations,  while  Curve  B  predicts  a  late  onset  with  low  severity  followed  by 
low  oscillations.  Curve  C  is  middle  grounds  between  B  and  D. 

Case  2 

The  imperfection  of  case  1  is  eliminated,  i.e.  y=0  in  Eqn.  24.  Using  a  fixed  mesh 
of  20x20  elements  with  vo=10  m/s,  At=2.0xl0‘9,  t2/At=6(XX),  and  solving  the  case  of  a 
softening  solid  (EjsIO^Eq),  we  obtain  curve  A  in  Fig.  12.  Although  there  is  no 

imperfection  in  this  case,  a  region  of  high  effective  strain  starts  developing  at  -u/ho=0.10 
and  acquires  the  shape  shown  in  Fig.  13a  for  the  last  calculated  point,  •u/ho=0.12.  The 
mode  of  localization  is  quite  different  than  that  of  case  1  (Fig.  9a)  and  occurs  at  a  later 
stage  of  loading:  -u/ho=0.10  as  compared  to  0.06  for  case  1.  Thus,  in  this  problem,  an 
imperfection  is  not  needed  to  achieve  localization  and  the  pattern  of  shear  bands  is  quite 
different  than  with  an  imperfection.  The  plot  for  the  SP  error  at  -u/ho=0.12  (Fig.  13b) 
indicates  high  error  in  the  localized  region.  This  will  lead  to  refinement  within  the  whole 
shear  band. 

The  problem  was  also  solved  starting  with  an  initial  mesh  of  20x20  elements,  and 
using  the  SP  error  measure  to  guide  the  adaptivity  with  e^  =  0.02,  e^'^=  0,  At=1.0xl0"9, 
ti/At=12(X)0,  t2=ti,  a=l/15,  and  hinin=  ho/80.  The  response  curve  B  is  shown  in  Fig. 
12  and  the  corresponding  adaptive  mesh  is  shown  in  Fig.  14. 

Curves  A  and  B  in  Fig.  12  are  the  response  curves  for  the  same  problem,  the  only 
difference  being  in  the  size  of  the  smallest  elements  in  the  shear  band.  The  curves  are 
identical  up  to  the  onset  of  localization  after  which  Curve  B  shows  a  steeper  stress  drop 
followed  by  higher  oscillations;  this  result  is  similar  to  case  1. 


Case  3 

In  this  case  we  consider  a  softening  solid  with  e|=200eQ  and  impose  a  velocity 
Vo=30  m/s.  This  case  is  identical  to  the  one  analyzed  by  Belytschko  et  al  (1992)  using  a 
fine  mesh  with  128x128  elements.  We  start  with  an  initial  coarse  mesh  (16x16  elements) 
and  use  the  SP  error  measure  to  guide  the  adaptivity.  Fig.  ISb  shows  the  curve  of  the 
average  stress  versus  end  displacement  corresponding  to  the  adaptive  mesh  shown  in  Fig. 
15a.  This  curve  is  in  good  agreement  with  the  one  for  the  fine  mesh  (128x128) 
represented  by  discrete  points  in  Fig.  15b.  The  adaptive  mesh  suggests  the  existence  of  a 
shear  band.  This  can  be  verified  by  plotting  the  profile  of  the  effective  strains  across  the 
band  along  the  normal  direction  at  -u/ho=0.12  (Fig.  15c).  This  profile  is  also  in  good 
agreement  with  the  results  for  the  fine  mesh. 

2.6.  CONCLUSIONS 

The  following  error  estimates  have  been  studied  in  one-dimensional  localization 
problems:  the  Zienkiewicz-Zhu  global  L2  stress.  Babushka-Rheinboldt  residual  criterion, 
the  interpolation  criterion  of  Diaz  et  al  (1983),  and  the  strain-projection.  The  last  is  a  new 
error  criterion  proposed  herein  which  is  based  on  an  L2  projection  of  strains. 

It  was  shown  that  error  indicators  such  as  the  momentum  equation  residual 
(developed  by  Babushka  &  Rheinboldt  (1978)  for  linear  elastic  statics)  and  the  L2 
projection  on  stresses  (developed  by  Zienkiewicz  &  Zhu  (1987)  for  linear  elastic  statics) 
are  not  effective  for  localization  problems.  They  fail  to  indicate  the  need  for  refinement  in 
the  area  of  localization,  where  the  strain  gradients  are  very  large  but  the  stress  is  nearly 
uniform. 

An  interesting  finding  from  the  adaptive  two  dimensional  calculations  is  that  in  the 
dynamic  problem,  imperfections  interact  with  the  amplification  due  to  wave  propagation 
effects  to  set  the  shear  band  pattern.  In  static  problems  with  homogeneous  stress  fields, 
the  imperfection  is  critical  in  setting  the  localization  and  the  length  scale  (i.e.  the  internal 
structure)  of  the  shear  band.  Thus  if  a  large  imperfection  is  placed  at  the  center,  a  pattern 
of  shear  bands,  with  one  passing  through  the  middle  develops.  Otherwise,  a  different 
pattern  of  shear  bands  develops. 

It  is  also  of  interest  to  note  that  while  nonuniform  meshes  in  wave  propagation 
problems  lead  to  additional  errors  due  to  spurious  reflection  at  interfaces  between  regions 
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with  different  element  sizes  (Holmes  and  Belytschko  (1976)),  in  adaptive  methods, 
varying  element  sizes  do  not  detract  from  accuracy.  The  reason  for  this  is  evidently  the 
tendency  of  the  fine  mesh  to  move  with  the  wave  front,  so  that  the  wave  front  does  not 
pass  through  regions  of  different  element  sizes.  In  fact,  the  coarser  elements  away  from 
the  wave  front  are  beneficial  because  they  filter  out  spurious  high  frequency  components. 

The  results  obtained  here  make  us  quite  optimistic  about  the  potential  of  h- 
adaptivity  in  explicit,  transient  finite  element  programs.  The  implementation  of  h- 
adaptivity  in  this  class  of  programs  is  not  difficult,  since  slave  nodes  which  arise  at  mesh 
gradations  are  easily  handled  in  an  explicit  setting.  Furthermore,  the  overhead  associated 
with  adaptivity  is  not  large,  and  a  more  accurate  solution  can  be  obtained  with  the  same 
computational  resources.  The  regions  of  localization  cannot  be  known  prior  to  a  run,  and 
the  use  of  a  uniform  mesh  with  the  refinement  necessary  to  resolve  shear  bands  accurately 
requires  excessive  computational  resources.  For  general  purpose  programs,  the  addition 
of  h-adaptivity  is  more  difficult,  but  nevertheless  the  advantages  are  promising. 
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Captions 


Fig.  2.1.  (a)  Fission-fusion  cycle  for  a  four-node  quadrilateral  and  a  two-node  bar,  (b) 
neighbor  configuration  for  a  four-node  quadrilateral  and  a  two-node  bar,  (c)  transition 
from  mesh  a  to  mesh  b  by  fission  of  element  2. 

Fig.  2.2.  Repeat  algorithm  for  an  initial  mesh  with  four  elements. 

Fig.  2.3.  Stress-strain  response  at  various  constant  strain  rates  (A,B,C)  and  a  variable 
strain  rate  (D). 

Fig.  2.4.  Comparison  of  (a)  strain  and  (b)  stress  profiles  for  a  fixed  80  element  mesh 
(80),  with  adaptive  solutions  using  Zienkiewicz-Zhu  error  criterion  (ZZ)  and  strain- 
projection  error  criterion  (SP). 

Fig.  2.5.  Error  profiles  for  one  dimensional  localization  problem  using:  (a)  strain 
projection  (SP),  (b)  interpolation  (IN),  (c)  Zienkiewicz-Zhu  (ZZ),  and  (d)  Babushka- 
Rheinboldt  (BR). 

Fig,  2.6.  Evolution  of  mesh  with  time  guided  by  (a)  Zienkiewicz-Zhu  and  (b)  strain 
projection  error  measures. 

Fig.  2.7.  Block  subjected  to  compression. 

Fig.  2.8.  Average  stress  versus  end  displacement  curves  for  the  block  using  a  hardening 
material  (Curve  A)  and  a  softening  material  (Curves  B,  C,D). 

Fig.  2.9.  Effective  strain  surfaces  at  -u/ho=0.12  corresponding  to  (a)  Curve  C  and  (b) 
Curve  A  in  Fig.  2.8. 

Fig.  2,10.  Strain  projection  error  surfaces  at  -u/ho=0.12  corresponding  to  (a)  Curve  C 
and  (b)  Curve  A  in  Fig.  2.8. 

Fig.  2.1 1.  Adaptive  mesh  evolution  corresponding  to  Curve  D  in  Fig.  2.8. 


Fig.  2.12.  Average  stress  versus  end  displacement  curves  for  a  softening  material  with  no 
imperfection. 

Fig.  2.13.  (a)  Effective  strain  and  (b)  Strain  projection  error  surfaces  corresponding  to 
Curve  A  (Fig.  2.12)  at  -u/ho=0.12. 

Fig.  2.14.  Adaptive  mesh  corresponding  to  Curve  B  (Fig.  2.12.). 

Fig.  2.15  (a),  adaptive  mesh,  (b).  Average  stress  versus  end  displacement  curve  for  a 
softening  solid  and  (c).  Effective  strain  across  the  shear  band  along  the  normal  direction. 
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Fig.  2.1.  (a)  Fissicm-fusion  cycle  for  a  four-node  quadrilateral  and  a  two-node  bar 
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(b)  neighbor  ccmfiguration  for  a  four-node  quadrilateral  and  a  two-node  bar 
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(c)  transition  from  mesh  a  to  mesh  b  by  frssion  of  element  2. 
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Fig.  2.2.  Repeat  algorithm  for  an  initial  mesh  with  four  elements. 
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Fig.  23.  Stress-strain  response  at  various  constant  strain  rates  (A,B,Q  and  a  variable  strain  rate  (D). 
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Fig.  2.4.  Comparison  of  (a)  strain  and  (b)  stress  profiles  for  a  fixed  80  element  mesh  (80),  with  adaptive  solutions  using 
Zienkiewicz-Zhu  error  criterion  (ZZ)  and  strain-projection  error  criterion  (SP> 
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Fig.  2.5.  Error  profiles  for  one  dimensional  localizauon  problem  using;  (a)  strain  projection  (SP).  (b)  interpolation  (IN), 
(c)  2Uenkiewicz-Zhu  (22),  and  (d)  Babuska-Rbeinboldt  (BR) 
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